function [T_synth failures]= rs_decode(R,n,k,b)
%
% T_SYNTH = RS_DECODE(R,n,k,b) 
%
% Reed Solomon decoding using the Euclidean Algorithm. 
% It works for both when (n-k) is odd and even. 
%
% R  needs to be a GF object to be decoded.
% n  is the number of symbols in the encoded codeword. This number can be
%    less than 2^m-1 in case of a shortened RS code.
% k  is the number of symbols in the decoded message. 
% b  is the factor in (X-alpha^B)*(X-alpha^(B+1))*...*(X-alpha^(B+N-K-1)),
% defined in the generator polynomial, where b is an integer. If it is not
% provided, b is assumed to be 1 by default.
%
% This function uses rs_init function in order the initialize the lookup
% tables that are saved as global variables  and rs_syndromes function to
% calculate the syndrome values of the given received polynomial. 
%
% Tansu Demirbilek
% ademirbi@mc.com
% Mercury Computer Systems

%% These values are calculated from the input variables. 
m = R.m;
N = 2^m-1; 
f = N-n;            % Code shortening factor
K = k+f;    % (N,K,m) for the original (unshortened) RS code
N_Tries = size(R,1);

% Calculate t depending on whether (n-k) is even or odd
if (mod(n-k,2)==0)
    t = (n-k)/2;
else
    t = (n-k-1)/2;
end

if(nargin == 4)
    % Use a default of b=1 like the Matlab does
    b=1;
end
p = R.prim_poly;

%% Check/Create the global GF lookup tables
gf_name = sprintf('GF_%d_%d_%d_%d_%d',N,K,m,b,p);
eval(sprintf('global %s',gf_name));

if(isempty( eval(gf_name)))
    rs_init(N,K,m,p,b);
end

% Create the local copies of the global variables
eval(sprintf('g = %s.g;',gf_name));
eval(sprintf('alpha_to = %s.alpha_to;',gf_name));
eval(sprintf('alpha_inv = %s.alpha_inv;',gf_name));
eval(sprintf('log_alpha = %s.log_alpha;',gf_name));

% Calculate the Syndrome vector
N_Tries = size(R,1);
S = rs_syndromes(R,n,k,b);
R = [zeros(N_Tries,f) R];

%% Calculate the Error Locator and Magnitude Polynomials
% (using Euclidean Algorithm)
% For euclidean algorithm formulation, check out:
% http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
% http://www.dsprelated.com/showmessage/22571/1.php


% Calculate the Error Locator and Magnitude Polynomials

% Initialize L and W matricies
% We are assuming that v<=t where v is the number of errors
L = gf(zeros(N_Tries,2*t+1), m,p);
W = gf(zeros(N_Tries,t), m,p);
list = []; % List of erroneous codewords
failures = []; % List of received codewords with decoding failures. 

for j=1:N_Tries
    % Only decode if syndrome vector is non-zero
    if (~isequal(S(j,:),gf(zeros(1,n-k), m,p)))
        % Keep a list of received codewords that was received in error
        list = [list j];
        
        % Allocate memory for Q and R matricies considering the worst case.
        % Each row of Q and R corresponds to a quotient and remainder calculated at
        % each iteration step.
        Q = gf(zeros(2*t+1, 2*t+1), m,p);  % Quotients Matrix
        D = gf(zeros(2*t+1, 2*t+1), m,p);  %
        A = gf(zeros(2*t+1, 2*t+1), m,p);  % Remainders Matrix

        d = 2*t; % Just an initial value so that the while loop starts
        % Set the initial values
        A(1,:) = gf([1 zeros(1, 2*t)], m,p); % Set the Dividend to be x^2t initially
        A(2,:) = polyclean(S(j,:), 2*t+1);
        D(1,:) = gf(zeros(1, 2*t+1), m,p); % [0 ... 0] i.e., d0 = 0
        D(2,:) = gf([zeros(1, 2*t) 1], m,p); % [0 ... 0 1] i.e., d1 = 1
        i=0;
        while (d >=t) % Keep iterating until the degree of \Omega 
                      % (Error magnitude polynomial) is less than t
            i = i+1;
            % The line below is the polynomial division and the most
            % cycle-intensive part of the algorithm
            
            % Divide A_i by A_(i+1)
            [Q_temp,R_temp] = deconv(polyclean(A(i,:),0),polyclean(A(i+1,:),0));
            % Add zeros to the higher order values so that R_temp will fit
            % into the matrix. Also get the degree of the remainder, to see
            % if we can stop the algorithm or not 
            [R_temp, d] = polyclean(R_temp,2*t+1);
            % Append the remainder into the A matrix
            A(i+2,:) = R_temp;
            % Keep track of Quotients 
            Q(i,:) = polyclean(Q_temp,2*t+1);
            
            % There might be an overflow at this point. If that happens,
            % Matlab will produce an error
            % The last D polynomial will be the 
            % Error Locator Polynomial L(x)
            D(i+2,:) = polyclean(conv(Q(i,:),D(i+1,:)),2*t+1) + D(i,:);
        end

        % Scaling factor is D_j(0), get the actual values of L(x) and W(x) by
        % dividing by this scaling factor
        L_deg = polydegree(D(i+2,:));
        if (D(i+2,2*t+1) ~= 0)
            if( (L_deg<=t+1) && (L_deg>=1) && (polydegree(A(i+2,:))<=t) )
                L(j,:) = polyclean(D(i+2,:)./D(i+2,2*t+1),2*t+1); % Error locator polynomial
                W(j,:) = polyclean(A(i+2,:)./D(i+2,2*t+1),t); % Error magnitude polynomial
            else % Decoding Failure
                failures = [failures j]; % Mark the code as decoding failure
                % i.e., there are more than t errors in received symbol,
                % cannot decode
                L(j,:) = gf(zeros(1,2*t+1),m,p);
                W(j,:) = gf(zeros(1,t),m,p);
            end
        else
            if( (L_deg<=t+1) && (L_deg>=1) && (polydegree(A(i+2,:))<=t) )
                L(j,:) = polyclean(D(i+2,:),2*t+1); % Error locator polynomial
                W(j,:) = polyclean(A(i+2,:),t); % Error magnitude polynomial
                % L(j,:), W(j,:)
            else % Decoding Failure
                failures = [failures j]; % Mark the code as decoding failure
                % i.e., there are more than t errors in received symbol,
                % cannot decode
                L(j,:) = gf(zeros(1,2*t+1),m,p);
                W(j,:) = gf(zeros(1,t),m,p);
            end
        end
    end

end

% Correct the Errors

% Assign the T_synth to R. Error vectors will be added to it. Codewords
% with no error and decoding failures will remain untouched. 
T_synth = R;

% T_synth = gf(zeros(N_Tries,n),m,p);
for j=1:length(list) % only process those that have non-zero syndrome (have error)
    ind = list(j);

    if(isempty(find(failures==ind))) % Only attempt to correct the error 
        % if the received codeword was NOT marked as a decoding FAILURE

        %%% Calculate Error Locations %%%
        
        % Solve the Error Locator Polynomial (by CHIEN SEARCH)
        % To find the error locator polynomial roots, basically try all of the
        % alpha values one by one. This is a trivial search since one needs to try
        % `n` number of alpha values (from alpha^0 ... alpha^(n-1))

        % L is evaluated at [alpha^(-1) ... alpha^(-14) alpha^(-15)]
        % Zeros give the roots of the L (error locations)
        % For N>n (Shortened codes), if an ErrorLocation is greater than n it
        % should be interpreted as a decoding error.
        ErrorLocations = find(polyval(L(ind,:),alpha_inv) == 0);

        %%% Solve the Error Magnitudes %%%
        % First, obtain L'(X_j^(-1)) = L_1 + L_3*X_j^(-2) + L_5*X_j^(-4) + ...
        L_prime = L(ind,1:end-1);
        L_prime(end-1:-2:1) = 0; % set even terms to zero

        %%% Calculate Error Magnitudes %%%
        if (polyval(L_prime,alpha_to(mymod(ErrorLocations*-1,N))) ~=0)
            Y = alpha_to(mymod(ErrorLocations*(1-b),N))...
                .*polyval(W(ind,:),alpha_to(mymod(ErrorLocations*-1,N)))...
                ./polyval(L_prime,alpha_to(mymod(ErrorLocations*-1,N)));

            E_synth = gf(zeros(1,n), m,p);
            % E_synth is indexed as (n-ErrorLocations) because the polynomials are set
            % as: [x^n-1 x^n-2 ... x^1 x^0]
            E_synth(mymod(n-ErrorLocations,N)) = Y;

            % If the code is shortened, only generate error symbols for the
            % first n symbols. For future reference: if there's an error symbol
            % calculated by the decoder at a location larger than n, it means
            % that the decoded symbol is in error and needs to be discarded.
            if(length(E_synth) > n)
                failures = [failures ind];
                %E_synth = E_synth(end-n+1:end);
            else
                %%% Correct the Received vector %%%
                T_synth(ind,f+1:end) = R(ind,f+1:end) + E_synth;
            end

        else % Undefined division by zero. Cannot decode...
            failures = [failures ind];
        end
    end

end

T_synth = T_synth(:,f+1:end);
failures = sort(failures);
